Superfluid and Mott Insulating shells of bosons 
in harmonically confined optical lattices 
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Weakly interacting atomic or molecular bosons in quantum degenerate regime and trapped in har- 
monically confined optical lattices, exhibit a wedding cake structure consisting of insulating (Mott) 
shells. It is shown that superfluid regions emerge between Mott shells as a result of fluctuations 
due to finite hopping. It is found that the order parameter equation in the superfluid regions is 
not of the Gross-Pitaeviskii type except near the insulator to superfluid boundaries. The excitation 
spectra in the Mott and superfluid regions are obtained, and it is shown that the superfluid shells 
posses low energy sound modes with spatially dependent sound velocity described by a local index of 
refraction directly related to the local superfluid density. Lastly, the Berezinskii-Kosterlitz-Thouless 
transition and vortex-antivortex pairs are discussed in thin (wide) superfluid shells (rings) limited 
by three (two) dimensional Mott regions. 

PACS numbers; 03.75.Hh, 03.75.Kk, 03.75 Lm 
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I. INTRODUCTION 

The recent experimental discovery of Bose-Mott insu- 
lating phases in optical lattices has generated an explo- 
sion of research in the ultra-cold atom community (see [ij 
for a recent review), and has helped to merge two major 
branches of physics: atomic-molecular-optical and con- 
densed matter physics. Most experiments thus far have 
relied on measuring the momentum distribution of the 
atoms after switching off the trap confining the atoms 
to infer the existence of a superfluid to insulator transi- 
tion [2, 3,, i4]. However, very recently, two experimental 
groups 0,01 have used spatially selective microwave spec- 
troscopy to probe in situ the superfluid-to- insulator tran- 
sition of ^^Rb in a three dimensional (3D) optical lattice 
with a harmonic envelope. In these experiments, the shell 
structure of the Bose-Mott insulating states was revealed 
for very deep lattices. Regions of filling fraction n = 1 
through n = b {n = 1 through n = 3) were mapped in 
the MIT |5| (Mainz |6|) experiment. Their observations 
in three-dimensional optical lattices lead to the confirma- 
tion of the Mott-insulating shell structure consisting of 
"Mott plateaus" with abrupt transitions between any two 
sucessive shells as proposed in two dimensional (2D) [7] 
optical lattices. 

One of the next frontiers for ultra-cold bosons in op- 
tical lattices is the search for superfluid regions separat- 
ing Mott-insulating shells. Eventhough the Mott shell 
structure was determined recently using microwave spec- 
troscopy 0, Q , any evidence of superfluid shells has re- 
mained elusive. Thus, in anticipation of the next exper- 
imental breakthrough, we study 2D and 3D optical lat- 
tices of atomic or molecular bosons in harmonically con- 
fining potentials, and show that between the Mott regions 
of filling fraction n and n + 1, superfluid shells emerge 
as a result of fluctuations due to finite hopping, and ex- 
tend our previous work on this topic ^ . This finite hop- 



ping breaks the local energy degeneracy of neighboring 
Mott-shells, determines the size of the superfluid regions 
and is responsible for the low energy (sound) and vortex 
excitations. In addition, we find that the order param- 
eter equation is not in general of the Gross-Pitaeviskii 
type. Furthermore, in 3D optical lattices, when super- 
fluid regions are thin (nearly 2D) spherical or ellipsoidal 
shells, we obtain bound vortex-antivortex excitations be- 
low the Berezinski-Kosterlitz-Thouless (BKT) transition 
temperature which is different for each superfluid 

region. Finally, we propose the use of Laguerre-Gaussian 
and Bragg spectroscopy techniques for the detection of 
superfluid shells. 



The remainder of the manuscript is organized as fol- 
lows. In Sec. nil we present the Bose-Hubbard Hamil- 
tonian in a harmonic trap. In Sec. IIIIl we analyze the 
emergence of an alternating superfluid and Mott shell 
structure, comparing two different approaches involving 
non-degenerate and nearly degenerate perturbation the- 
ory. We also obtain the order parameter equation, the 
characteristic sizes of the superfluid and Mott regions, 
the local filling fraction and the local compressibility. We 
find that the superfluid order parameter that emerges 
between two Mott shells is not of the Gross-Pitaeviskii 
type, except very close to the insulating boundaries. In 
Sec. IIV[ we describe quasiparticle and quasihole excita- 
tions in the Mott regions and quasiparticle, sound and 
vortex excitations in the superfiuid regions. In Sec. |Vl 
we present a possible experiment using Gauss-Laguerre 
beams and Bragg spectroscopy, which can be performed 
in order to visualize the existence of superfluid shells. 
Lastly, we state our conclusions in Sec. |VT1 
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FIG. 1: (Color online) Phase diagram of the homogeneous 
Bose Hubbard model. For the inhomogeneous Bose Hubbard 
system, the dashed (red) line indicates the values of the local 
chemical potential fir that the system exhibits from /ir — 2.5(7 
(FiglU or fj,r — 2.2U (FigO at the center of the trap to /ir = 
at the edge of the trap for t — 1.25 x 10~^U. It indicates 
how the system exhibits an alternating stucture of Mott and 
superfluid shells 



II. BOSE-HUBBARD HAMILTONIAN 

To describe the physics of alternating insulating and 
superfluid shells of atomic or molecular bosons in optical 
lattices, we use the lattice Bose-Hubbard Hamiltonian 
with a harmonic potential described by 



H 



-t ^ 4cr+a + 
r,a 



u 



V/Xr4cr (1) 



where fir — fJ- — V{r) is the local chemical potential, 

V{r) = np{p/af/2 + n,{z/a)y2 (2) 

is the harmonically confining potential, a is the lattice 
spacing, and cj is the creation operator for a boson at 
site r. Here, the vector r = {x,y,z), and the distance 
p = \J + y^. The anisotropic confining potential V{v) 
becomes isotropic when J7p = fi^. Furthermore, t is the 
hopping parameter and JJ is the interaction strength, 
which we assume to be repulsive (positive). The sum in 
the first term on the right hand side of Eq.[T]is restricted 
to nearest neighbors. The harmonic trap makes the 
system inhomogenous and introduces interesting proper- 
ties which are absent in homogenous Bose-Hubbard sys- 
tems For the inhomogeneous systems, it is useful 
to define a local chemical potential [i-r — [i—V(y) which 
plays an essential role in the emergence of shell struc- 
tures. 

For bosons in a homogeneous optical lattice, the Bose- 
Hubbard model without the harmonic confining potential 
can be used. In this case, it is well known that an inte- 
ger number of particles on each site or r + a describes an 
insulating state for sufficiently large U {U ^ t). This is 
because the on-site interaction U makes it energetically 



unfavorable for a particle to move from one site to an- 
other. In this situation the^stem is in what is known as 
the Mott insulator phase [13|. However, for non-integer 
number of particle , the extra bosons can move more eas- 
ily, at a small energy cost, because its interaction energy 
is essentially the same on every site. For this reason, a 
system with non-integer number of bosons on each site 
is a superfluid at zero temperature Recently, sev- 

eral Mott-insulator shells of bosons were detected in op- 
tical lattices Q, which are inevitably inhomogenous, 
and thus require a model that includes the effects of the 
harmonic part of the confining potential. This can be 
modeled by the inhomogenous Bose-Hubbard hamilto- 
nian described in Eq. (jT]) . In this case, one can either be 
in a regime where the entire system is superfluid or in 
a regime where the system exhibits a shell structure of 
alternating Mott-insulator and superfluid regions. The 
existence of this shell structure has been shown numeri- 
cally 01 , but the lack of analytical progress has hindered 
a true understanding of the emergence and properties of 
these shells. 

A simple argument for the emergence of the shell struc- 
ture can be made by inspection of the standard phase 
diagram of the homogenous Bose-Hubbard model shown 
in Fig. [T] upon the substitution of /i — + /ir. In an inho- 
mogenous system, the effective local chemical potential 
/ir is spatially varying. Thus, in the regions where the 
number density n(r) is not an integer, a superfluid shell 
emerges, and in the regions where n(r) is an integer, a 
Mott-insulator shell appears. Eventhough, a simple argu- 
ment for the existence of the shell structure can be made, 
many questions need to be seriously addressed. For in- 
stance, what are the characteristic order parameter, di- 
mensions and excitations of each superfluid shell? Thus, 
we begin our presentation by discussing next the emer- 
gence of the Mott-insulator and superfluid shell struc- 
ture. 



III. EMERGENCE OF THE SHELL 
STRUCTURE 

A standard approach to analyse such bosonic systems 
is to use the Bogoliubov mean field approximation. How- 
ever, as shown in [12|, this approximation fails to predict 
the expected phase transition since it treats the inter- 
actions only approximately. Hence, instead of using the 
Bogoliubov approximation, we generalize a method found 
in the literature pTl [T^ . by introducing a local mean 
field theory that treats the interactions exactly and ap- 
proximates the kinetic energy of the atoms in the op- 
tical lattice. We introduce the local superfluid order 
parameter = (cr)- We can now construct a consis- 
tent local mean field theory by susbstituting the operator 

cJCr+a -> (4>Cr+a + 4(Cr+a) " (cj ) (Cr+a) , leading tO an 



effective local Hamiltonian 

Hf = i/o,n (r) - i V (Cr V;+a + 



■V'rV'r+a), (3) 
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which is diagonal in the site index r with 

Ho,n{r) = ^hriflr - 1) - /ir^r, (4) 

where — cjcr is the number operator. For t ^ 0, the 
sheh structure for Mott-insulating phases is revealed by 
fixing fir = n, to obtain the local energy 

-E'o,n(i") = -y"!"- - 1) - Mrn, (5) 

when (n — 1)U < /ir < nU . Since i?o,n+i(r) — £'o,n(r) = 
nil— fir, the change from a Mott shell with filling fraction 
n to n + I occurs at the degeneracy condition /ir = nU , 
which for a spherically symmetric potential happens at 
the radius 

Rc,n = a^Jrin/^l, (6) 

where fin — 2{fi~nU). The relation /ir = nU determines 
the shape and size of the boundary between the n and 
n + 1 shells. For instance, in the case of the anisotropic 
potential of Eq. [2] the same condition leads to ellipsoidal 
shells 

with principal axes Op = a^JVln/^p, and = a^J^n/^z- 
However, near this region of degeneracy, fluctuations due 
to hopping introduce superfluid shells, as discussed next. 

A. Continuum approximation 

In this section, we show how the hopping term in 
Eq. (121) affects the ground state energy of the system. 
Qualitatively one can see that the kinetic energy con- 
trolled by t lifts the degeneracy of the system at /ir = nU 
and in the process introduces a superfluid order param- 
eter in a region of finite width depending on parameters 
n, and U . 

To obtain analytical insight into the emergence of su- 
perfluid shells, we make first a continuum approximation 
through the Taylor expansion 

V'(r + a) = ^{r) + aj9j?/'(r) + ^aia-,c'i9-,V(i"), (8) 

where repeated indices indicate summation. This should 
be true in the limit where the order parameter is 
smoothly varying at a length scale much greater than 
the unit lattice spacing a. Under this approximation the 
effective local Hamiltonian becomes 

Hf - i/o,„(r) - CrA*(r) - 4A(r) + A(r) (9) 

where A(r) = ztip{r) + ia^V^V(r) reflects the amplitude 
for the creation of a single boson excitation at position 
r, while the last term 

A(r) = i[A(r)V'*(r) + A*(r)V(r)] (10) 



reflects the local mean-field energy shift. Here z — 2d 
is the number of nearest-neighbor sites (coordination 
number) depending on the lattice dimension d. Within 
this approximation, a simple analytical treatment of the 
emergence of superfiuid shells is possible, as discussed 
next. 



B. Nearly Degenerate Perturbation Theory 

We focus our attention now on the Mott regions with 
integer boson filling n and n + 1 and the superfluid shell 
that emerges between them. In the limit where U ^ t 
we can restrict our Hilbert space to number-basis states 
\n) and |n -I- 1) at each site. Any contribution of other 
states to the local energy will be of the order of t'^/U. The 
hopping term in Eq. ^ affects the ground state energy of 
the system by removing the local degeneracy of i?o,n+i (r) 
and £'o,Tt(r) at /ir — nil. To illustrate this point, we use 
the continuum approximation described above and write 
the Hamiltonian in Eq. [9] in the matrix form 

^eff _ f ^o, „(r)+ A(r) -V7rTTA(r) \ 
^[ -V^rTTA*(r) i?o,„+i(r) -I- A(r) J ' 

where A(r) and A(r) is defined in Eq. [TOl Notice that 
t 7^ has two effects. First, it changes the local ener- 
gies i?o.n(r) and -Eo,n+i(r) of the Mott shells n and n+1 
through A(r). Second, it mixes the two Mott regions 
through the off-diagonal term \/n + lA(r) and its her- 
mitian conjugate. Thus, the physics near the boundary 
between the n and n+1 Mott regions is described by an 
effective local two- level system with diagonal (A(r)) and 
off-diagonal {y/nl^A{r)) perturbations. 
The eigenvalues of Eq. pT|) are given by, 

E±ir) = E,ir) ± v'(i?,(r))2 + (n + 1) |A(r)^ (12) 

where Es{r) = [£'o,„+i(r) + Eo,n{r)] /2 + A(r) is propor- 
tional to the sum of the diagonal terms, and Ed{r) = 
[£^o,ri+i(r) - Eo^ri{r)] /2 = {nil - /ir)/2 is proportional 
to their difference. These local eigenvalues are illustrated 
schematically in Fig. [2l where the local energies £'o,n(r) 
and £'o.n+i(r) are shown together with the eigenvalues 
E±{r). The radii Rn.± indicated in the figure correspond 
to the locations where A(r) = 0. 

Notice that -^-(r) is the lowest local energy, leading 
to the ground state energy E — J drE-{r), where L 
is the characteristic linear dimension of the system. 

The order parameter equation (OPE) is determined by 
minimization of E with respect to '0*(r) leading to 

-(-) (n+lM. + «-V-)A(r) . (^3) 
2^\E,{rW + {n + l)\A{r)\' 

Notice that the OPE is not of the Gross-Pitaeviskii (GP) 
type, since the superfiuid regions emerge from local fiuc- 
tuations between neighboring Mott shells. The zeroth 
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FIG. 2: (Color online) Schematic plot of local energies Eo,n{r) 
and i5o,n+i(r) showing the degenerate radius Rc,n- This local 
energy degeneracy is lifted by the the presence of a finite 
hopping t, which leads to an avoided level crossing shown 
as dashed dark-grey (red) curve, and to the emergence of 
a superfiuid region with inner radius R„.- and outer radius 

Rn,+ 



order solution of this equation with A(r) = ztip{r) leads 
to the spatially dependent order parameter 



+ 1 {nU-iir) 



(14) 



Since |'0(r)| > 0, it implies that \nU — fir\ < {n+l)zt in 
the superfluid region. Thus the inner radius Rn,- and the 
outer radius Rn.+ of the superfluid shell between the n 
and n + 1 Mott regions is obtained by setting |'0(r)|^ = 
leading to 



Rn,± — Rc 



2zt{n+l) 



'1± 



i?2 ■ 



(15) 



for a spherically symmetric harmonic potential V{r) = 
r2(r/a)^/2, where Rc,n is defined above. Our notation to 
describe the Mott-superfluid boundaries is illustrated in 
Fig.H 

Equation [15] shows explicitly that t splits the spatial 
degeneracy of the n and n-\-l Mott shells at r = i?c,n (or 
/ir = nU) by introducing a superfluid region of thickness 



A_R„ — R„ 



Rn.-- 



(16) 



In the case of non-spherical harmonic potential V{r) = 
Q.p{p/af' /2 + Vtz{z/ay' /2 the shell regions are ellipsoidal 
instead of spherical. Notice that Ai?„ depends strongly 
on filling fraction n, the ratio zt/rt and the chemical po- 
tential ji through i?c,ji- 

In addition, the local filling fraction 



ri(r) = — 



nU ^ 



2zt{n + l) 



(17) 



FIG. 3: (Color online) Schematic plot of superfluid regions, 
shown in gray (red), and the Mott regions, shown in white. 
The radius Rc,n separates Mott shells with filling factor n and 
n-\-l for t = Q. The superfluid regions have inner radius Rn.- 
and outer radius Rn,+ and emerge between Mott shells n and 
n + 1 for non-zero t. 



in the same region interpolates between n+1 for r ^ Rn.- 
and n for r ^ Rn.+i while the chemical potential /i is fixed 
by the total number of particles N = J drn{r). 

In Fig. H n(r), |^/'(r)r, and i?„,± are shown for the 
Mott and superfluid regions, for t — 1.25 x IQ^'^U, and 

= 6 X 10~^U, and fj, = 2.5U. For these parameters, 
three Mott and three superfluid shells emerge. It is very 
important to emphasize that in the superfluid regions 
the order parameter |'0(r)P is not identical to n{r) since 
the OPE equation is not of the Gross-Pitaeviskii type. 
While the order parameter ^/;(r) vanishes at the bound- 
aries Rn,± between the superfiuid and Mott shells, and 
reaches the maximum value |V'(r)|ma2; = (" + when 
fir = nil , the average particle density n(r) interpolates 
harmonically between Mott shells n+1 and n, having 
the average value of n + 1/2 when /ir — nU . Further- 
more, we show in Fig. [5] that when t — 1.25 x 10~^J7, 
and = 6 X 10~®C/, and n = 2.2U the central shell is su- 
perfluid, and the order parameter has a local minimum at 
the origin of the harmonic trap. As fi increases, this min- 
imum is reduced to zero at a critical value, and the n = 3 
Mott shell emerges. This property of the emergence of a 
Mott shell from the center of the trap by supressing the 
superfluid order parameter at that region is a generic fea- 
ture, and can be inferred directly from the phase diagram 
of Fig. [T]via the substitution fi ^ fir- 

The local compressibility 



K(r) 



1 



9n(r) 
dfi ~ 2zt{n+ 1) 



(18) 



of the superfluid shells is non-zero, in contrast to the in- 
compressible (k = 0) n and n+1 Mott shells for r < i?„._ 
and r > Rn.+ , respectively. The local compressibility in- 
dicates the presence of large (small) excitation energies 
in the Mott (superfluid) regions, as discussed later. 
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FIG. 4: (Color online) a) Shell structure of Mott and su- 
perfluids regions in a 2D square optical lattice with harmonic 
envelope as a function of radius r/a for t ^ 0. The superfluid 
regions are shown in red (gray) whereas the Mott regions are 
shown in white. The black circles indicate the Mott bound- 
aries Rc,n at t = 0. b) The local filling factor n(r) is shown in 
solid lines for t ^ and in dashed lined for {t = 0). The red 
curve (solid gray) shows the local superfluid order parameter 
\Tp{r)f. The parameters are Q = 6xlO~^U, t = 1.25x 10"^;7 
and n = 2.5(7. 



FIG. 5: (Color online) a) Shell structure of Mott and super- 
fluids regions in a 2D square optical lattice with harmonic 
envelope as a function of radius r/a for t ^ 0. The superfluid 
regions are shown in red (gray) whereas the Mott regions are 
shown in white. The black circles indicate the Mott bound- 
aries Rc,n at t = 0. b) The local fllling factor n(r) is shown in 
solid lines for t ^ and in dashed lined for {t — 0). The red 
curve (solid gray) shows the local superfluid order parameter 
|V'(r)p. The parameters are Q = 6xlO~^U, t = 1.25x 10"^!7 
and fi = 2.2(7. 



We note that only near the edges of the superfluid re- 
gion (where r « Rn.± and ?A(r) « 0) a direct expansion of 
the OPE (Eq. [T5|) leads to the effective Gross-Pitaeviskii 
equation 



1 



2moff 



+ l/eff(r) + .geff = 0. (19) 



Here, h = 1, rrieff = l/2a^t is exactly the boson band 
mass due to the optical lattice, T4ff(r) = \nU — fir\/{n + 
1) — zt is the effective potential, and gctf — 2zt/{n + 
1) is the effective interaction. Notice that T4ft(r) < 
and vanishes at the boundaries Rn,± of the superfluid 
region since \nU — /Zr| = zt(n + 1) there. Furthermore, 
f^cff = zt/{n + 1) is small in comparison to U, indicating 
that the superfluid near the edges is weakly interacting, 
and more so as the Mott index n increases. When t 
(moff oo) then |'(/'(r)|^ = —Vcs/gcS leading to the 
correct limiting behavior of Eq. [T3]near r « Rn,±- 

As discussed above, our nearly degenerate perturba- 
tion theory method provides a good description of the 



emergence of superfluid shells. This method is a gener- 
alization of the perturbation theory developed in [111 Il2l | 
for the uniform case, which was believed not to be ex- 
tendable to describe the emergence of superfluid in har- 
monically confined optical lattices [l4|]. Our method, de- 
scribed above and in Q, which relies on direct diago- 
nalization of an effective two level system, thus provides 
the connection between the mean-field pseudo-spin pic- 
ture described in , and the perturbative approach de- 
scribed in [l^ . 

Next, to clarify when the standard theory developed 
in [ill . Il^ breaks down, we compare the results from 
our nearly-degenerate perturbation theory with the non- 
degenerate case. 



C. Non-degenerate perturbation theory 

When the local energies for Mott phase n with energy 
£'0,n(r) — Un{n— \)/2 — ^^n, and Mott phase n-\-\ with 



6 



energy i?o^„+i(r) = U{n + l)n/2 — ^r{n + 1), are away 
for the degeneracy region — nU, then the correction 
to the local energy i?o,n is 



EL'' - E 



\{n\V\m)\^' 
Eq.. 



En 



(20) 



where |m) denotes the unperturbed wave function of 
the excited state with eigenvalue i?o.m- Here, V — 
— CrA*(r) — cjA(r) couples only to states with one more 
or one less atom than in the ground-state, and represents 
the perturbation to the Hamiltonian ffo,n(r) defined in 
Eq.m The, fourth order correction can also be calculated 
using higher-order perturbation theory , and leads to 
the Ginzburg-Landau energy 



E^^ao + A(r) + a2|A(r)p + a4|A(r)|^ 



(21) 



The coefficients ao, a2 and 04 are all functions of param- 
eters n,U,^r, and t. The term oq corresponds to the 
unperturbed energy i?o,n(r) described in Eq. [5] and as- 
sociated with the unperturbed Hamiltonian Ho^n defined 
in Eq. 01 while A(r) is the energy shift shown in Eqs. [5] 
and 1101 which contains spatial derivatives of the order 
parameter field V'(r). 

The second order coefficient 



n+1 



1 



0-2 



U{n—l) — fir iir — Un zt' 



(22) 



determines the existence of superfluid regions, when 02 < 
0, while the fourth order coefficient 



04 



n{n — 1) 



[U{n - 1) - fir] [Ui2n - 3) - 2/Xr] 



(nr + l)(n + 2) 



[fi-UnY[2fir-U{2{n+l))] 
( n n+1 



\J(n—V) — \ir Hr — Un^ 
n n+1 



{U{n - 1) - firf (Air - Unf 



(23) 



is always positive and is essentially identical to the ho- 
mogeneous non-degenerate limit [IJ] except for the re- 
placement of fj, —^ fj,r- However, the inherent inhomo- 
geneity of the trap potential manifest itself in the local 
energy for the superfluid regions (Eq. [^T|) through the 
spatial dependence of the coefficients 00,02, and a„ and 
through the spatial derivatives of ip{r) contained in A(r). 
Thus, we define the local Ginzbug-Landau energy dif- 
ference AEn = En — ao, which, in terms of the order 
parameter tp, becomes 

AE„ = -ta^il;*V^il} + a2zH'^ |V(r)|^ -I- a^zH'^ y)\'^. 

(24) 

Here, we retained only terms up fourth order in -0, and 
up to second order in derivatives of -0. 



Minimizing AEn to zeroth order in the spatial deriva- 
tives of ij}, leads to 



|V(r)r 



02 



2a4z2i2 ' 



(25) 



and setting |'0(r)| = (or 02 = 0) leads to the local 
chemical potential 



2/i± = U{2n -l)-zt± y/U^ -2U{2n+l)zt + zH^, 

(26) 

which determines the inner and outer radii of the Mott 
shell with filling n, where the order parameter ip{r) = 
vanishes. Solving Eq. [26] gives the smaller radius of the 
Mott shell with filling n, 



Rn.+ — Rc.n\ / 1 



and the larger radius 



2zt{n + 1) a2 



'1- 



2ztn 



d2 



(27) 



(28) 



calculated to order t'^/U. Here, Rn.+ > Rc.n, and 
Rn-i.- < Rc,n-i, where Rc.n radius of the Mott shell 
with filling n when f = 0, as defined in Eq. [6l Notice 
that the size of the Mott region ARn,Mott = Rn-i- — 
Rn.+ always decreases with increasing t from its value 
Rc.,n~i — Rc,n at t = 0, showing that superfluid regions 
emerge at the expense of shrinking the Mott insulator 
shells. 

The corresponding radii defining the superfluid region 
between the Mott shells with filling factors n and n + 1 
are precisely Rn,+ and Rn,-, where Rn,± are defined in 
Eq. [TS] Thus, to order P/U'^, the thickness of the super- 
fluid shells ARn is again given by Eq. [161 and we recover 
the results obtained from our degenerate perturbation 
theory. 

Furthermore, minimization of AEn (Eq. I24p with re- 
spect to ijj reduces to the same Gross-Pitaeviskii equation 
described in Eq. [THjuear the boundaries Rn,+ and Rn,-, 
where the order parameter V'(r) vanishes. The mapping 
near the boundaries is a^t l/2TOoff, a2z'^t^ Vcsi'^), 
and a/^z^t'^ — > Coff- 
in Fig. [Hi we compare the results of the non-degenerate 
perturbation theory with that of our nearly degenerate 
perturbation theory for small values of t/U. The pa- 
rameters used are same as in Fig. [Jj and we show the 
results only for the innermost superfluid shell, between 
n — 2 and n = 3 Mott shells. Notice that the non- 
degenerate perturbation theory is correct only very close 
to the boundaries between the superfluid and Mott re- 
gions where iP{t) is small, but fails to describe the su- 
perfluid phase at the center of the superfluid shell corre- 
sponding to the degeneracy region. 

Having, discussed the non-degenerate perturbation ap- 
proach and its breakdown, we analyze next the excitation 
spectrum in the Mott and superfluids regions. 
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FIG. 6: (Color online) The squared amplitude of the super- 
fluid order parameter |i/)(r)p is shown as a solid line (red) 
for the nearly degenerate case, and as a dashed line (gray) 
for the non-degenerate case. Notice that near the boundaries 
where i/'(r) « 0, both methods agree and describe the su- 
perfluid accurately. However, at the center of the superfluid 
shell, the non-degenerate method breaks down. The parame- 
ters are same as in Fig. |4l and we show the order parameter 
for the inner most superfluid shell, between n — 2 and n = 3 
Mott shells. 



IV. EXCITATIONS IN MOTT AND 
SUPERFLUID REGIONS 

In this section, we discuss relevant excitations in the 
Mott and superfluid regions, and we use the method of 
functional integrals to obtain quasiparticle and quasihole 
excitations in the Mott shells and sound and vortex ex- 
citations in the superfluid shells. 



A. Quasiparticle and quasihole excitations 
in the Mott regions 

The excitation spectrum in the Mott shells can be ob- 
tained using the functional integration method [l^ for 
the action {n = ks = 1, P = l/T) 



5[ct, c] = / dr V [4 ,a,Cr,r + H] 

Ja 



(29) 



leading to the partition function Z = 
/PctPcexp (-S'[ct,c]). In each Mott shell we in- 
troduce a Hubbard- Stratonovich field ^E* to take into 
account fluctuations due to the presence of finite hop- 
ping, and integrate out the bosons (c^ , c) leading to an 
effective action 



2Ci;,kk' 



to quadratic order in ^'^ and where k, k' are momen- 
tum labels, Lu are Matsubara frequencies, and 



l + fik 



iuj — El (r) iuj ~ E2 (r) 



with Ei{r) — nil — /ir, E2{r) — {n — 1)U — fir, and 
The poles of Gkk' r) are found upon the analytical 
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FIG. 7: (Color online) a) Quasiparticle Eqp (solid line), quasi- 
hole Eqh (dashed line), and Mott gap Eg (red) energies for 
k = versus r/a. b) Sound velocity for the outermost super- 
fluid ring versus r/a. Same parameters as in Fig. |4] 



continuation iu — lu + i5, leading to the local excitation 
energies 



^^1 

2 2 



- (4n + 2)ekC/ + f/2 



where the +(— ) sign labels quasiparticle (quasihole) ex- 
citations. The energy to add a quasiparticle is Eqp = w^, 
while the energy to add a quasihole is Eqh — —lo^. 

Figure [7^ shows the quasiparticle and quasihole en- 
ergies of the Mott phase as a function of position r. 
The energy cost to create a quasiparticle (quasihole) 
is minimum (maximum) at the trap center and in- 
creases (decreases) radially, while the Mott-Hubbard gap 
Eg = mink(-Egp + Eq}-^) is large and independent of r. 
Thus, Eg = mink\/ek - 4(n + 2)ekt^ + C/^, which leads 
to the final result Eg = \/{ztY - 4(n + 2)ztU + U^. 
This expression indicates that it is easier to create a 
quasiparticle-quasihole excitation inside higher n Mott 
shells. The horizontal solid (red) line in Fig. [7^ also shows 
this tendency. Since this gap is reduced with increasing 
n the Mott-insulator gets weaker and thus more suscep- 
tible to superfluid fluctuations. Notice that Eg vanishes 
when zt/U ~ (2n + 1) ± 2^n{n -I- 1), but the physical 
solution for the critical value of zt/U corresponds to the 
largest value, and thus the plus (-f ) sign. 

The large excitation energies in the Mott regions give 
away to low energy excitations in the superfluid regions 
as discussed next. 
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B. Excitations in superfluid regions 

While the creation of quasiparticle-quasihole excita- 
tions in the Mott regions can be energetically costly, the 
creation of single quasiparticle excitations, sound waves 
and vortex excitations in the superfluid shells are more 
easily accessible. 

Single quasiparticle excitations: The single quasiparti- 
cle excitation energy can be read off from Eq. [T^ since 
the first excited state of is i?+(r), and the ground 
state is i?_(r). The local energy difference is A£' = 
i?_i_(r) — -E'-(r) is independant of r. This can be seen 

from the expression A£' = 2\J\Ed{Y)\'^ + {n + l) |A(r)|^ 
when the approximate result A(r) « tz^j{v) is used in 
combination with Eq.[T4land with the definition Ed{v) = 
[£'o,„+i(r) — iJo,n(r)] /2 — {nU — fXr)/'^- The final answer 
is AE sa {n + l)tz, which indicates that energy cost of 
adding a quasiparticle in the superfluid shell is small in 
comparison to the cost of adding a quasiparticle in the 
Mott region. 

It is also interesting to analyze the eigenvectors of the 
local Hamiltonian defined in Eq. [121 The eigenvector 
corresponding to -E+ is 



-Ea+AE/2 
Vn+lA*(r) 
1 



(30) 



and reduces to the vector (0, 1)"^ corresponding to the 
energy En of the Mott phase with filling n, when A ^ 0. 
The eigenvector corresponding to E- is 



\E) = 



1 



(31) 



-E^-AE/2 



and reduces to the vector (1,0)"^" corresponding to the 
energy En+i of the Mott phase with filling n 4- 1, when 
A^O. 

However, the most interesting excitations in the super- 
fluid regions are collective in nature. We will discuss next 
the collective sound excitations and later the appearance 
of vortices. 

Sound velocity: The excitation spectrum of collective 
modes in the superfluid region can also be calculated 
using the functional integral method. First we intro- 
duce the Hubbard-Stratonovich field ip which now cor- 
responds to the order parameter in the superfluid re- 
gion. Second we use an amplitude-phase representation 
■0(r, t) = |'0(r, r)| exp [iip{r, t)] and apply the nearly de- 
generate perturbation theory described earlier to inte- 
grate out the boson fields and c. Thus, we obtain the 
phase-only effective action 



SeS 



1 

21? 



drdr [k(9i-'p)^ + pijdi(fdjip] (32) 



to quadratic order in the phase variable for the super- 
fiuid region between the n and n + 1 Mott shells. Here, 



we assumed that [-(/"(r, t)| is r-independent at the sad- 
dle point. The coefficient k is the compressibility of the 
superfluid described in Eq. [181 aiid the local superfluid 
density tensor 



{n+l)ta^ F{\il^\,n, t) 



2ta^\ij\^S,j, (33) 



2 GM,n,t) 
where the numerator of the first term is 
F{\^\,n, t) =. t {Az\ib\^d,j + 4|^|VVl% - 2d,\^P\dj 
and the denominator of the first term is 



Gm,n,t) - ^{nU~fir)y^+{n + lWTm, 

with the function 

Tm = z^\^\' + 2\^PM^\ + {W^\^\f. 

The complex structure of the superfluid density tensor is 
a direct consequence of the non-Gross-Pitaeviskii nature 
of the order parameter equation (Eq. I13p . 

Insight can be gained into the structure of the local 
superfluid density tensor by neglecting the gradient terms 
involving 1-01, which in combination with Eq. [141 produces 
a local superfluid density tensor 



Ps{r) = pu = 2ta^\i){Y) 



(34) 



which vanishes at the Mott boundaries Rn,±- This local 
superfluid density tensor has been described previously 
in Refs. [?] and however the more general expression 
shown in Eq. [331 goes beyond the mean fleld approxima- 
tion presented in the pseudo-spin description [1^ . In the 
present approximation, the resulting wave equation has 
the form 



= 0, 



(35) 



leading to a local sound velocity c(r) = ^JpJ^r)/K, 
which in terms of the order parameter reads c(r) — 
2-^/(n + l)zta|^/'(r)|. The local speed of sound has it max- 
imal value at |'0(r)|maa; = (%/?! + l)/2, and the superfluid 
region behaves as a medium of continuous index of refrac- 
tion 



X(r) 



^max '\/tI -\- "T 

^^2\My)\- 



(36) 



Notice that x(r) ^ oo at the Mott boundaries where 
|V'(r)| — 0, indicating that the sound waves of the super- 
fluid do not propagate into the Mott regions. A plot of 
the local sound velocity is shown in Fig. [Tti for the su- 
perfluid region between the n = 1 and n = Mott shells. 
From the phase-only effective action for the superfluid re- 
gion, we can also investigate the vortex excitations, which 
are discussed next. 

Vortices and antivortices: Next, we explore vortex so- 
lutions in two cases where spontaneous vortex-antivortex 
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pairs can appear as indicators of the Berezinski- 
Kosterlitz-Thouless (BKT) transition [§, • Case I cor- 
responds to a 3D system, where the superfluid regions 
are very thin Ai?„ <^ Rc,n, leading to essentially a two- 
dimensional superfluid in curved space. Case II corre- 
sponds to a 2D system, where the superfluid regions are 
thick rings Ai?„ ~ Rc,n, leading to essentially a two- 
dimensional superfluid subject to boundary conditions. 

In a flat space two-dimensional system stationary vor- 
tex solutions must satisfy f \/ip ■ dl = 27rm, where 
TO = ±I,±2, ... is the vorticity (topological charge) and 
V(/3 is the superfluid velocity. The standard vortex so- 
lution in cylindrical coordinates is Vip = mO/r, and the 
corresponding free energy per unit volume is 

T^^^ j dvp,{v){V^)\ (37) 

This situation is analogous to a two-dimensional linear 
dielectric material where the displacement field is 

D X z = e(r)E, (38) 

with dielectric function e(r) — I/[27rp(r)]. Notice that 
the dieletric function diverges at the superfluid bound- 
aries, since Ps(r) — ^ in those regions. In this language 
T is identical to the electrostatic energy per unit volume 

C^el = ^/^rD.E. (39) 

In general the solutions for several vortices (antivortices) 
can be obtained from 

VA V</5 = 27rz^m,(5(r-r,), (40) 

i 

where is the location of the vortex (or antivortcx) of 
vorticity rrii. 

In case I the superfluid state appears below Tbkt ~ 
7rps(r = Rc.n)/2, where ps = Psjc?' has dimensions of 
energy. In this limit ps(r = Rc,n) = (n -|- l)i/2 and the 
critical temperature of the superfluid shell between the n 
and n + \ Mott regions depends on the index n. Notice, 
however that such estimate is reasonable only when the 
radius of the shells are sufficiently large. The solution for 
a vortex-antivortex (VA) pair in curved two-dimensional 
space is 

^(. = i?,„, ^, 0) = arctan i^^^f"^^^) 

where i?c,n = a-\/2(/i — nU)/il is the radius of the su- 
perfluid shell, and h is the VA size. A three-dimensional 
view of the velocity field V</3(0, 0) is shown in Fig. [S] 
When the superfluid shell has a small thickness Ai?„ then 
Tbkt ~ nt{n -\- I)Ai?„/6a, while the vortex-antivortex 
pair has an approximate solution of the same form as 
above which interpolates between ip{r = Rn_-,9,(j)) and 
ip{r = Rn,+,9,(j)). 



In case II the superfluid regions are rings bounded by 
Rn.- and Rn,+ , and one can use the Coulomb gas anal- 
ogy described above, conformal mapping techniques and 
proper boundary conditions to obtain vortex-antivortex 
solutions. The creation of vortex-antivortex pairs are en- 
ergetically quite costly when A_R„ <C Rc,n, due to strong 
conflnement effects of the boundaries, thus we do not 
expect a BKT-type superfluid transition to occur until 
Ai?„ is substantially large (~ Rc.n)- Only in this limit, 
we expect a BKT transition with Tbkt ~ '^(Ps) /2, where 
{ps)/2 is the surface area average of Ps(r). 




FIG. 8: Three-dimensional view of a vortex-antivortex pair 
in 2D superfluid shell separating two 3D Mott regions. 

Having analyized the excitation spectra in the insula- 
tor and superfluid regions, we discuss next some possible 
experimental detection schemes for the superfluid shells. 

V. DETECTION OF SUPERFLUID SHELLS 

Evcnthough the use of tomographic microwave tech- 
niques has allowed the detection of several Mott re- 
gions @, Q, it has been quite challenging to detect su- 
perfluid shells. In this section, we propose a possible 
experiment for the detection of superfluid shells through 
the use of Gauss-Laguerre and Gaussian beams followed 
by Bragg spectroscopy. The idea is that Gauss-Laguerre 
and Gaussian beams can tranfer angular momentum to 
the atoms in superfluid phase without transfering lin- 
ear momentum, and that Bragg spectroscopy can detect 
their existence since the technique is only sensitive to 
the velocity of the atoms in the superfluid phase. This is 
because the Mott insulator regions do not absorb the an- 
gular momentum because of the presence of a large gap 
in the excitation spectrum. Next, when the frequencies 
of the beams used in Bragg spectroscopy are correctly 
chosen some of the rotating atoms in the superfluid re- 
gions acquire extra momentum and are kicked out of their 
shells and can be imaged. 

The Gauss-Laguerre technique has been sucessfully 
used to rotate superfluid sodium atoms ^"^Na in a mex- 
ican hat potential [l^, [l^. The technique used was a 
two step process, where initially a Gauss-Laguerre beam 
transfered both linear and angular momentum to atoms. 



and a second Gaussian beam canceled the linear momen- 
tum transfer leaving each trapped atom with exactly one 
unit of angular momentum. In this manner, the atoms 
participating in superfluidity rotate without dissipation, 
and when the trap is released, the superfluid region ex- 
pands, but does not fill the center of the cloud, thus main- 
taining a toroidal structure throughout the time of flight. 
Furthermore, Bragg spectroscopy was used to detect the 
sense of rotation of vortices within a vortex lattice of 
sodium BEG (TTj and to observe the persistent flow of 
Bose-condensed Sodium atoms in toroidal traps [l^ . 

The Gauss-Laguerre technique and Bragg spectroscopy 
may also be used to detect superfluid shells of bosons in 
harmonically confined optical lattices. However, the sit- 
uation here is slightly different because of the existence 
of Mott shells and multiple superfiuid regions. To illus- 
trate this we discuss the simpler case of a nearly two- 
dimensional configuration, where the harmonic trap is 
very tight along the z-direction and loose along the x- and 
y- directions and only two superfluid shells are present. 
Upon application of the Gauss-Laguerre technique along 
the z- direction, angular momentum transfer occurs es- 
sentially to the atoms in the superfluid phase, imposing a 
rotating superfluid current with a well defined superfluid 
velocity profile, while the Mott regions remain unchanged 
due to their large gap in the excitation spectrum. The 
angular momentum transfer is chosen to occur along the 
z direction generating the superfluid currents in the x- 
y plane, and the amount of angular momentum trans- 
fer is assumed to be Ti for each atom that absorbed a 
Gauss-Laguerre photon with ^ = 1. Next we use Bragg 
spectroscopy with counter-propagating beams along the 
x-direction to detect the sense of rotation and determine 
the regions of superfluid rings with well defined velocities. 
If there is sufficient optical resolution and signal-to-noise 
ratio, then the experiment may be performed in situ, 
otherwise the Bragg spectroscopic measurements can be 
performed in time-of-fiight. 

A schematic plot of the detection of superfluid shells 
using Bragg spectroscopy can be found in Fig. [9l which 
shows two superfluid shells rotating counter-clockwise 
and two counter propagating beams applied along the 
x direction. As can be seen in Fig. [5]the right(left)-going 
beam has frequency uj{uj') and linear momentum k(k'). 
The bose atoms undergo a transition from the internal 
state with energy to the internal state e/. In the fol- 
lowing analysis, we retain h, instead of setting h = 1. 
Applying momentum conservation, we can easily obtain 
that the final linear momentum of the atoms in the su- 
perfluid region is p/ = pi -I- ?i(fc -I- /e')x, in terms of the 
initial linear momentum of the photons hkx. and -hk'i.. 
Thus, the Bragg beams transfer a net linear momentum 
h{k + fc')x to the atoms which satisfies the energy con- 
servation condition 

h(LU'-uj') = £f-e,-vMk + k') + —- '-!- (41) 

2m 

In the equation above Vx is the component of the veloc- 
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FIG. 9: (Color online) Schematic plot for the detection of su- 
perfluid shells using Bragg spectroscopy. The angles 6i and 62 
indicate the locations of strongest momentum transfer from 
the Bragg beams (large green arrows) to the rotating super- 
fluid shells of radii Tii and R2- The gray arrows indicate the 
sense of rotation of the superfluid shells. 



ity Vi = Pi/m along the x direction. Notice that Vx can 
also be written as Vx = Vi sin 9 where 9 is the angle be- 
tween the Bragg beams and the velocity of the atoms 
in the superfiuid shell and Vi = Pi/m is the magnitude 
of Vi. For an atom carrying one unit of angular momen- 
tum, the superfluid velocity is Vj = h9/mr. Therefore, 
within a superfluid shell at position r ^ R atoms get a 
linear momentum kick of h(k + fc')x when the velocity 
Vx — h sin 9 /mR satisfies the condition given in Eg. 1411 
This leads to two Bragg angles 9 = — sin"^ (mRvx/h), 
and TT — 9 for each superfluid shell. As can be seen in 
Fig. [5] the Bragg angles are 9i and n — 9i for the outer 
superfluid shell labelled by and are 6*2 and n — 92 
for the inner superfluid shell labelled by i?2- The re- 
gions with same velocity Vx are identified by the condition 
sin6'i/i?i = sin6'2/i?2 = mvx/h. Once these atoms are 
kicked out of their respective superfluid shells, they form 
two small cloud, which can be detected by direct imaging. 
As mentioned before, if there is sufficient optical resolu- 
tion and signal-to-noise ratio, then the experiment may 
be performed in situ, otherwise the Bragg spectroscopic 
measurements can be performed in time-of-flight. In fact 
these clouds also carry the information of the sense of ro- 
tation of the superfluid shells and their velocity profile, 
and one can in principal extract that information from 
the images. 

Having discussed our proposal for the experimental de- 
tection of superfluid shells in harmonically confined op- 
tical lattices, we present next our conclusions. 
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VI. CONCLUSIONS 

We studied 2D and 3D optical lattices of atomic or 
molecular bosons in harmonically confining potentials, 
and showed that between the Mott regions of filling frac- 
tion n and n + 1, superfluid shells emerge as a result 
of fluctuations due to finite hopping. We found that 
the presence of finite hopping breaks the local energy 
degeneracy of neighboring Mott-shells, determines the 
size of the superfluid regions as shown in Fig [51 and 
is responsible for low energy (sound) and vortex exci- 
tations. In addition, we obtained an order parameter 
equation which is not in general of the Gross-Pitaeviskii 
type, except near the boundaries separating the super- 



fluid from the insulating regions. Furthermore, we ob- 
tained bound vortex-antivortex solutions (as shown in 
Fig M below the Berezinski-Kosterlitz-Thouless (BKT) 
transition when superfluid regions are thin (nearly 2D) 
spherical (or ellipsoidal) shells. Finally, we discussed 
that the emergence of these superfluid regions should be 
detectable using a combination of Gauss-Laguerre and 
Bragg spectroscopy techniques. 
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